Inversion of spinning sound fields 



Michael Carley, 
Department of Mechanical Engineering, 
University of Bath, 
Bath BA2 7AY, 
England; 

Electronic mail: m.j.carley@bath.ac.uk 
November 3, 2008 



Abstract 

A method is presented for the reconstruction of rotating monopole source distribu- 
tions using acoustic pressures measured on a sidehne parallel to the source axis. The 
method requires no a priori assumptions about the source other than that its strength 
at the frequency of interest vary sinusoidally in azimuth on the source disc so that the 
radiated acoustic field is composed of a single circumferential mode. When multiple az- 
imuthal modes are present, the acoustic field can be decomposed into azimuthal modes 
and the method applied to each mode in sequence. 

The method proceeds in two stages, first finding an intermediate line source derived 
from the source distribution and then inverting this line source to find the radial vari- 
ation of source strength. A far-field form of the radiation integrals is derived, showing 
that the far field pressure is a band-limited Fourier transform of the line source, es- 
tablishing a limit on the quality of source reconstruction which can be achieved using 
far-field measurements. The method is applied to simulated data representing wind- 
tunnel testing of a ducted rotor system (tip Mach number 0.74) and to control of noise 
from an automotive cooling fan (tip Mach number 0.14), studies which have appeared 
in the literature of source identification. 

1 Introduction 

This paper describes a method for determining rotating source distributions from acoustic 
measurements. This is a problem which has been examined by a number of researchers, with 
many [1-6] considering the problem of estimating the amplitudes of the acoustic modes at 
the termination of a circular duct, as in the case of aircraft engines. The motivation for these 
studies has usually been to determine the source terms in their own right, in order to find 
the source mechanisms responsible for the noise or to improve noise control measures, but a 



second application has been in developing models which can be used to predict the acoustic 
field. This prediction model can be used in active noise control [7, 8] or in using near- 
field measurements taken in a wind-tunnel to make far-field predictions of noise radiated 
by aircraft in flight [1,9]. This gives rise to two different, though related, problems: the 
first is the determination, to within some tolerance, of the acoustic source; the second is 
the determination of the acoustic source to within a tolerance sufficient to give accurate 
predictions of the acoustic field at points other than the measurement positions. 

This paper considers a model problem for the recovery of a rotating source distribution 
from a set of measurements along a sideline, a line parallel to the source axis. The question 
of how to position microphones, and how many to use, features in the analysis of many 
researchers. Typical microphone configurations have included 3 microphones at 120 angu- 
lar positions [1], 91 microphones on a fixed polar array [6], 18 microphones rotating over 
20 positions [2] and 21 microphones located on a fixed arc [4], depending on the experimen- 
tal facilities used and the fidelity of results required. Recent work on engine noise has also 
included the use of sensor arrays mounted inside or on the engine. Examples are the use 
of 100 pressure sensors on the surface of the intake[10] and simulations of an array of 150 
microphones mounted on the wall of an engine duct [11]. In these cases, the methods used are 
described as 'beamforming' and come from the class of techniques used for source location 
rather than for source characterization. 

In this paper, we present an inversion technique which uses data from a linear arrange- 
ment of microphones to recover the details of a distribution of monopoles on a disc. This 
corresponds to the problem of thickness noise of a propeller or other rotor [12], sound from a 
baffled circular piston [13] or to soimd radiated by the termination of a circular duct, when 
the Rayleigh approximation is valid [14, 15]. The only assumption made is that the source, 
and the acoustic fleld, have a known sinusoidal variation in azimuth — no assumption is re- 
quired about the form of the radial variation of the source nor is a far-fleld approximation 
needed. The resulting method is applied to simulated data using parameters characteristic 
of problems to which identification methods have been applied in the past. 

The source recovery technique which is developed here is based on the measurement 
techniques used in wind tunnel measurement of aerodynamic sources[l-9] but could also 
be viewed in the more general framework of sound source reconstruction in other areas of 
acoustics[16] and, in particular, in relation to cylindrical nearfield acoustical holography 
(NAH)[17, 18] where measurements are taken on a cylindrical surface surrounding a source 
region and then forward projected to find the acoustic field elsewhere in space, or back- 
projected to find the acoustic quantities which characterize the source. The method of this 
paper shares some similarities with NAH but differs in incorporating known information 
about the source geometry and azimuthal dependence. 

2 Inversion of spinning sound fields 

The acoustic fields to be considered in this paper can all be viewed as being generated by 
sources distributed over a disc. Figure 1 shows the arrangement of the problem. We use 
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Figure 1: Coordinate system for radiation prediction 



cylindrical coordinates (r, 6*, z) with an acoustic monopole source distributed over the disc, 
< r < 1, 2; = 0, non-dimensionalizing all lengths on disc radius. At a single frequency cu, 
the acoustic field p is given by [12]: 



p{r,9,z,u;) = s{n,9,)^d9,ndn, (1) 



"1 r^n ^]kR 

'0 Jo 

with wavenumber k = lu/c, c the speed of sound and source-observer distance: 

^ + - 2rri cos{e - Oi) + z^. 

The source term s(ri, Oi) can be decomposed into a series of azimuthal modes with sinusoidal 
variation: 

00 

n=— 00 

which, upon insertion into Equation 1 with the transformation 9 — 9i ^ 9i, yields: 

00 pi /•27r ^{kR-nOi) 

p{r,9,z,u;)^ J2 / ^n(ri) . „ d9,ndn, (2) 

^t^^ Jo Jo 47ri? 

with R being redefined: 

= + rf — 2rTi cos 9i + z'^. 

The acoustic field at a frequency uj is thus a sum of azimuthal modes, each of which is 
directly generated by a corresponding azimuthal mode on the source disc. The aim of the 
inversion procedure is to recover the source function(s) Sn(^i) given as input some acoustic 
pressures measured in the field. The nature of these measurements will depend on the type 
of source being studied. 

There are two main categories of problem which will be considered: rotating sources such 
as propellers and fans and ducted sources where the duct termination can be considered a 
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disc-shaped source. For a source rotating at angular frequency Jl, the radiated field contains 
only harmonics of frequency nQ. Furthermore, if the source strength is steady in the rotating 
reference frame, there is only one azimuthal mode, of order n, present at each of these 
frequencies. This means that the acoustic field of Equation 2 reduces to: 

/•I /•27r ^(fciJ-nei) 

p(r, e, z, nn) = ei"M / s^in) ^ ^ d^i n dn. 
Jo Jo 47ri{ 

The properties of the acoustic field are largely controlled by the rotor speed and, in 
particular, the tip Mach number which, for a source of unit radius, is Mi = fl/c. When the 
source rotates supersonically, Mt > 1, the acoustic field is dominated by the source around 
the sonic radius r* — 1/Mi [19]. When Mt < 1, the blade tip is the dominant region and, 
in the far-field, its radiation is exponentially stronger than that from inboard regions [20]. 
This means that the measured field is effectively the field radiated by the tip and recovering 
the details of the source at smaller radii will be difficult. On the other hand, if the aim is 
to accurately compute the acoustic field at a new set of points, it may well be sufficient to 
capture only the source at the tip. 

The structure of the rotating field has been studied using model solutions [21-23] and 
the role of the sonic radius has been clarified. The field is made up of a segmented near 
field which undergoes a transition around r*. In 'tunneling' across this transition region, the 
sound field decays exponentially, explaining the relatively weak field radiated by subsonically 
rotating sources. Supersonic sources have part of the source lying beyond r* so that they 
can radiate strongly into the field, without losing energy in tunneling through the transition. 
This transition region means that source recovery will always be a hard problem if only far- 
field data are available, a result which will be derived in §2.3 by considering the bandwidth 
of the spatial data in the far field. 

When the radiating system is a circular duct, the problem can be modeled by taking the 
noise source to be the duct termination. In that case, the source distribution is composed 
of the duct modes which have propagated to the end of the duct. The field inside a circular 
duct is composed of modes of the form Jnikmrj) exp[j (n6' — kzmnz)], with J„( ■ ) the Bessel 
function of the first kind, J'n{kmn) — and kzmn is an axial wavenumber [14]. When kzmn 
has an imaginary part, the mode decays exponentially in the duct and does not propagate 
to the termination. In any case, the source strength at the duct termination can be taken 
to be the acoustic velocity generated by the modes which do propagate and the radiated 
noise can be accurately computed over much of the field using a Rayleigh integral [3, 14] 
or a Kirchhoff integral over a wider range of polar angles [3,24]. The source can again be 
modeled as a circular disc with an azimuthally varying source term. A number of methods 
have been developed for the identification of the radiating modes [2-6] and have been found 
to be accurate and robust, considering the assumptions made in their development. 

The one extra difficulty in the duct case compared to the rotor problem is that the source 
at the duct termination may be composed of modes of more than one azimuthal order. In 
this case, there are procedures which use measurements at multiple angles to extract the 
modal amplitudes in the acoustic field. For example, a method has been presented which 
uses 360 measurements distributed over a semicircular 'hoop' to find the amplitudes of the 
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Figure 2: Arrangement of experimental measurements 



azimuthal modes radiated from a duct [2]. The hoop of microphones was then moved to find 
the modal amplitudes as a function of axial displacement z. 

From the known properties of rotating acoustic fields and established experimental tech- 
niques, it is clear that it is possible to measure and/or extract the complex amplitude of a 
single azimuthal mode radiated by a disc-shaped source. Indeed, if the source is tonal so 
that the modal content does not change with time, the measurements could, in principle, be 
performed with only two microphones, one fixed as a phase reference, and another moving 
along the sideline. 

2.1 Formulation 

Figure 2 shows the basic experimental arrangement. The input to the inversion method is the 
amplitude of a single azimuthal mode p(r, z) with r fixed. When the sound is generated by 
a steady rotating source, p(r, can be found by measuring the field on one sideline. When 
modes of different azimuthal order are present, the field must be measured on multiple 
sidelines of the same radius r, varying ^, and a decomposition procedure applied to find 
p(r, as discussed in the previous section. 

However the acoustic field may have been measured and processed, the sound radiated 
by one source mode of azimuthal order n at frequency u is found by integration over the 
source disc [12]: 



where the observer is positioned at (r, 0, z). The aim of the inversion algorithm is to recover 
the radial source distribution f{r\) from the field pressures p(r, z). 

To begin to recover f{ri), the first stage is to rewrite Equation 3 in a transformed coor- 
dinate system (r2, 02, z) centred on the measurement sideline. Figure 3. This transformation 
has been used in calculations of transient radiation from pistons [13, 25] and in studies of 
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Figure 3: Coordinate systems (ri, 6i, z) and (r2, 62, z). The (r2, O2, z) system is shown dashed 
and the thick line shows the region of integration over 62 in the transformed system. 



propeller noise fields [21-23, 26]. Transforming Equation 3 gives p(r, z) as an integral over a 
line source K{r,r2)'- 



P{'^^^)= I —s-K{f^^2)r2^r2, (4) 

Jr-\ ^ 

K{r,r2)^-j^^^^ e-J"^7(ri)d^2, (5) 
for observer positions with r > 1. The original coordinates (ri, 6*1) are related to (r2, O2) by: 



— r + r2 + 2rr2 cos ^2, (6a) 
_i r2sin^2 

9i = tan — . (6b) 

r + r2 cos 62 

so that the limits of integration in Equation 5 are given by setting ri = 1: 

= cos-^ (7) 
^ 2rr2 ^ ^ 

The function K{r, r2) depends only on the observer lateral displacement and is constant 
for all points on a sideline parallel to the source axis. The inversion method proposed is to 
measure p{r,z) at fixed r, invert Equation 4 to recover K{r,r2) and then use Equation 5 to 
recover /(ri). 



2.2 Inversion algorithm 

The first stage of the inversion procedure is to use measured sideline data to recover the 
source function K{r,r2). Noting the behavior of K at its endpoints, Equation 20, we write: 

K{r,r2) = [{r2 - {r - l)){r + 1 - r2)Y/^K'{r,r2). (8) 
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The integral of Equation 4 is discretized to give: 



1=1 ■' 



where: 

1/2 



Zj is the axial displacement of the jth measurement point and {t^^\ wf^"^) are the nodes and 
weights of an iV-point Gauss-Chebyshev quadrature rule of the second kind. 

Equation 9 can be written as a system of equations relating the vector of measured 
pressures p to the unknown vector of sources K': 

[A]K' = p, (10) 

A„~ir^r)yrK (11) 

In practice, the system will be over-determined, with the number of measured pressures M 
being greater than iV, the number of values of K' to be determined. At this stage, the 
system is solved for K', using some suitable method for ill-conditioned problems, with K 
being recovered from Equation 8. 

The second stage in determining the source distribution is to invert Equation 5 to recover 
/(ri). We proceed by approximating /(ri) as a sum of Legendre polynomials Pq(ri): 

Q 

f{n) = J2^Mri), (12) 

q=0 



SO that 



1 - /-Stt-^) 

^(^'^2) = irE^W(o) e-j"^^n(^i)d^2, (13) 

giving rise to the system of equations: 

[B]F = K, (14) 



^"=4^i« «-'"*n('-i)d9., (16) 

r2^r + tl ' . 

The integration is performed using a standard Gauss-Legendre quadrature. As before, this 
system can be solved using a method suitable for ill-conditioned problems and f{r\) recon- 
structed from the coefficient vector F. 
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2.3 Far- field limitations 



The integral of Equation 4 is identical to the exact integral of Equation 3. If we make the 
standard far-field approximations, we can establish some limit on the accuracy of reconstruc- 
tion possible using far field results. Expanding R to first order in r2: 



^ / \ T-, r 2 , 21 1/2 

n ~ no -|- . - - _ „ ^ 

so that: 



R^Ro + -^{r2-r), Ro= [r^ + z^] , 



p ^ / eJ^''''2/^°i^ (r, r2)r2 dr2, 

Rq Jr-l 



(16) 



which can be rewritten: 



P 



Rq 

a = kr/Ro, 



/oo 
^""-'Kir, r2)r2H{r2 - (r - l))H{r + 1 - r2) dr2, (17) 
-oo 



where H{ ■) is the Heaviside step function [27]. 

In the far field, Equation 17 shows that the measured pressure on a sideline is propor- 
tional to a band-limited Fourier transform of the source term K{r,r2). A reconstruction 
algorithm based on far field measurements can only recover components of K{r, r2) with 
spatial frequency < a < k, with components outside this frequency band being lost in 
tunneling across the transition region between the near and far fields. 

This result provides a link between XAH and the method of this paper. In XAH, a Fourier 
transform on the sideline data, i.e. Equation 17, is used to recover the coefficients of a field 
expansion in cylindrical wave-functions[17, 18].This leads to difficulties with finite aperture 
effects due to the periodicity enforced by the finite Fourier transform. In this algorithm, no 
use is made of the Fourier transform in the reconstruction procedure so that the shortcomings 
of the finite Fourier transform do not cause the spurious sources which appear in NAH. On 
the other hand, the discretization introduced by the finite number of samples on the sideline 
can lead to aliasing as in NAH and any other reconstruction procedure based on spatial 
sampling of the acoustic field. The implication of Equation 17 is that in order to use the 
information which is present on the sideline, the sampling rate must be such as to capture 
behavior up to wavenumber k. If the minimum sampling rate is taken to be twice per 
wavelength then the sideline measurements should be taken no more than ir/k apart. 



3 Results 

Two test cases have been simulated as a first test of the algorithm of §2.2. The first uses 
parameters characteristic of the CRISP ducted rotor tests[l] and the second models an 
automotive cooling fan which has been used in tests of noise control[7, 8]. In each case, 
the sound field p{z) is computed by integration of Equation 3. To simulate measurement 
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errors and background noise, a Gaussian random signal of amplitude emax|p| is added to 
the computed pressures before using them in the inversion scheme. 

For the ducted rotor test case, = 0.74, k = 7.4, n = 10, M = 128, N = 64, r = 1.125 
and < 2; < 4. The source /(ri) was synthesized by adding the first four duct modes of 
circumferential order n with a random phase so that the source was given by: 

4 

m=l 

where (pm is a random phase < 0^ < 27r. The number of measurement points N was 
chosen to be approximately equal to that in the CRISP tests where data were taken at 120 
points[l]. 

In the cooling fan case, M = 0.14, k = 0.83, n = Q, M = 16, N = 16, r = 1.25 and 
< z < 8. This time, the source used was /(n) = (1 - ri)^/^ as this gives reasonable 
physical behaviour near the blade tip[20]. Again, the number of sensor positions was chosen 
to be similar to that used in the original work: in this case, 17 microphones were used in the 
authors' source reconstruction experiments[7, 8]. 

The two test cases which have been chosen represent two realistic problems with quite 
different characteristics. The CRISP case is similar to many wind tunnel tests which aim to 
extract the acoustic source from in-field measurements: the source is quite high frequency 
and the tip Mach number is such that, although energy is lost in the transition to the far 
field, the acoustic field is quite strong and there is sufficient information to allow the source 
to be determined reasonably accurately. The low-speed cooling fan, however, presents a 
rather more difficult problem. Due to the low rotor speed, the field decays rapidly inside the 
sonic radius r* = 7.14 and the measured field has lost much of the content useful for source 
reconstruction. 

The inversion method has been implemented using Octave[28] and the Regularization 
Tools package of Hansen[29, 30]. Equation 10 is solved using Hansen's implementation of 
truncated singular value decomposition [31] with the regularization parameter automati- 
cally selected using the L-curve criterion [32]. The same technique is then used to solve 
Equation 14 and find f{ri). 

Two measures are used to assess the accuracy of the method. The first is to compare 
the recovered source g{ri) with the input f{ri). The second is to use g{ri) to compute the 
acoustic field q{r, z) at a new set of points and compare this field to p(r, z) computed using 
/(ri). This assesses the ability of the algorithm to 'project' measured data into the field. 

3.1 Source reconstruction 

The inversion algorithm has been run with zero added noise and with e = 10~^, equivalent 
to a maximum signal-to-noise ratio of 60dB. Figures 4-7 show the reconstructed source for 
the two test cases, with the source terms weighted on radius ri, as in the radiation integrals. 
The source reconstruction in the ducted fan case Figures 4 and 5 is quite good in both cases. 
With zero noise, it accurately reproduces the shape and amplitude of the input source. With 
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Figure 4: Ducted fan test case, source term, e = 0; solid and dashed lines, and 
circles and squares, ^{g) and '^{g). 
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Figure 5: Ducted fan test case, source term, e = 10 
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Figure 6: Cooling fan test case, source term, e = 0. 

added noise, the reconstruction is not quite as good, especially for inboard ri ^ 0.8 but the 
details of the source are captured quite well near ri — 1, the dominant region for radiation 
at this wavenumber. 

The cooling fan results, Figures 6 and 7 are not as good, probably because the number 
of sensors is quite small but also because the acoustic field is so much weaker than in the 
ducted fan case, due to the low rotor speed. As discussed in §2, sound from source regions 
inside the sonic radius decays exponentially as it radiates. Here the whole source lies inside 
the sonic radius, meaning that the acoustic field is composed largely of evanescent waves, 
making source reconstruction difficult. 

In Figure 6, the reconstructed source oscillates considerably at small radii, but the tip 
behaviour is very well captured. This might be expected: the tip is strongly dominant 
meaning that the recovery of the inboard source is very poorly conditioned. With noise 
added, the reconstructed source is smoother, although the amplitude is not found accurately. 
The tip behaviour, however, is again accurately computed. 

3.2 Field estimation 

Figures 8-12 compare the field computed using g{ri) to the real field p(r, z), near to (r = 2) 
and far from (r = 8) the source disc, for the e = 10^^ case. The results have been scaled 
on p(r, 0) to simplify comparison. Real and imaginary parts are shown separately as a check 
on the ability of the method to calculate the phase of the field, important in scattering 
calculations and in control. 

The ducted fan results, Figures 8 and 9, are very good. The phase has been accurately 
computed in the near and far fields and the amplitude error is about 10% of peak amplitude. 
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Figure 13: Cooling fan test case, reconstructed far-field noise, e = 0. 

or IdB. The directivity of the source is such that the field does not decay rapidly on the 
sideline, aiding the reconstruction technique. As a check that the method does converge to 
a correct result in the absence of noise, the reconstruction method has also been applied to 
data with e = 0. The recomputed far-field pressures are shown in Figure 10 and, as they 
should be, are very close to the correct data, indeed practically indistinguishable from them 
with the amplitude error at 2; = being 0.06dB. 

In the cooling fan case. Figures 11 and 12, the field decays rapidly and is reconstructed 
quite poorly. The shape and phase are roughly correct but the amplitude error is about 50% 
or 4dB. The error may be due to the form of the field or to the small number of sensors 
simulated. Note that although the amplitude of the reconstructed source is much less than 
that of /(ri), the reconstructed field amplitude is rather larger. This is due to the exponential 
dominance of the tip region as an acoustic source on subsonic rotors, mentioned in §2: the 
difference in tip gradient for g{r\) has made the computed acoustic field stronger than that 
found using 

In any case, given that the phase has been accurately computed, the result might still be 
useful in control applications where the phase of the control signal is important in cancelling 
the unwanted noise. Again, we present results with no added noise, Figure 13, and, here, 
the comparison is not as good as in the ducted rotor case. The shape of the field has been 
well captured but the amplitude is underestimated by about 3% or 0.26dB. 

3.3 Algorithm performance 

To assess the performance of the method when used with varying numbers of measurements, 
the source reconstruction method has been applied to the simulated CRISP data with M — 
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Figure 14: Error e versus number of sensor positions M and M/N: squares: M — N; circles 
M = 2N; triangles: M = AN. 




Figure 15: Condition number k of reconstruction matrices versus M, N = M: squares: [A] 
circles: [B]. 
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16,64,128,256 and M/N — 1,2,4. The error measure for the reconstructed source is the 
Loo norm: 

^ _ max|ri/(ri) - riff(ri)| 
°° max|ri(/(ri)| 

where the weighting with ri has been retained, corresponding to an area weighting of the 
error. The calculation has been performed with no added noise, to check the factors which 
contribute to the error in the reconstructed quantities. 

Figure 14 shows the variation of error with the number of sensors, as a function of the 
ratio M/N. The first obvious point is that, for this set of operating parameters, the error for 
M = 64 is very large when M/N = 4, while the method failed completely at M = 32. This 
appears to indicate that the source term cannot be well-approximated by only 16 terms in 
the expansion of Equation 12, an unsurprising result. 

More interesting is that for N/M = 1,2, the error decreases steadily as M increases, 
but then increases between M = 128 and M = 256. Figure 15 shows the condition number 
K, of the matrices [A] and [B] used in the inversion procedure, as a function of M, with 
N — M. As might be expected, the condition number of both increases with M, as the 
systems become more poorly conditioned. Machine precision on the computer used for the 
calculations is approximately 1/2^^. The condition number k(A) of the matrix used to 
estimate K{r,r2) is always greater than 2^° so that the first part of the inversion scheme is 
always ill-conditioned. The reason for the drop in accuracy past M = 128 seems to be the 
higher condition number of the second matrix used k(B). At M = 128, it rises above 2^^ 
and we conjecture that at this point the loss of precision in calculations is too great for the 
inversion method and the results begin to worsen. This is a function of the solver used and 
it may be that a different choice of regularization scheme would lead to better results. 

4 Conclusions 

A source reconstruction method for the inversion of spinning acoustic fields has been devel- 
oped and tested on two representative problems. It has been found that the method can 
work well, even with added noise, depending on the type of source to be identified. The 
method requires no a priori assumptions about the form of the source other than that it 
be circular and vary sinusoidally in azimuth. This makes it a useful intermediate between 
nearfield acoustical holography, where no information is assumed about the source except its 
approximate location, and other source identification methods which use assumptions about 
the location and spatial variation of the source to model its radiation characteristics. In 
the case of the method of this paper, it may be that when additional information about the 
source is available, such as its modal structure, users might be able to incorporate this in- 
formation into the technique to improve source reconstruction and/or to reduce the number 
of measurements required. 
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A End-point behaviour of K{r, r^) 



To establish the behavior of K{r,r2) near the endpoints of the integrand in Equation 4, we 
note that as r2 — >■ (r — 1), 6*2°'' — > tt and 6i — >• 0. When r2 — > (r + 1), 6*2°'* — >■ tt and 6i — >• tt. 
We examine the basic integral: 



K= — I ' e-j"^id^2. (18) 
47r ./fl(o) 



For 6*2°'* — TT and resulting small 9i: 

K^— ld^2. (19) 

Integrating, 

and using cos~^ x ^ {1 — x^)^/^ as x — > — 1 to give: 



,(o)^_ 2V2(^ + l-r2)V2(r2-(r-l))V2 



(2rr2''^'^^ 



2 

yields 



(.-1))./^ (20) 



with square root behaviour as r2 — (r — 1)+ and r2 — > (r + 1) . 
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